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Computing  Flow  through  Well 
Screens  Using  an  Embedded  Well 

of  Engineers®  Technique 

by  Hwai-Ping  Cheng 


PURPOSE:  The  purpose  of  this  Coastal  and  Hydraulics  Engineering  Technical  Note  (CHETN) 
is  to  document  a  computational  technique  developed  to  compute  the  flow  rates  through  the 
screens  of  groundwater  pumping  wells  or  relief  wells.  The  technique  was  developed  for  the 
three-dimensional  (3D)  groundwater  (GW)  flow  simulation  using  the  finite  element  (FE) 
method. 

BACKGROUND:  Accurate  estimation  of  flow  rates  through  well  screens  is  essential  in 
groundwater  modeling  and  may  impact  decisions  based  on  the  simulation  results  associated  with 
various  project  alternatives.  In  reality,  the  flow  rate  through  a  well  screen  can  vary  both  spatially 
and  temporally  due  to  heterogeneity  and  variability  present  in  the  surrounding  subsurface 
environment.  The  location  of  the  pump  within  a  pumping  well  may  also  impact  the  flow  rate 
distribution.  However,  a  uniform  distribution  of  flow  throughout  the  screen  length,  resulting 
from  the  total  pumping  rate  divided  by  the  screen  length,  has  been  often  assumed  and  employed 
to  characterize  groundwater  withdrawal  in  groundwater  modeling  due  to  its  simplicity.  This  may 
cause  inaccurate  model  calibration  and  lead  to  poor  decisions  as  a  result.  For  relief  wells  used  to 
protect  earthen  levees,  accurate  estimates  of  flow  rate  through  well  screens  under  various 
hydrologic  conditions  is  essential  to  the  design  of  both  the  relief  well  system  and  the  associated 
surface  conveyance  system  so  that  excess  groundwater  coming  out  of  relief  wells  shall  be 
diverted  efficiently,  effectively,  and  economically. 

In  this  CHETN,  a  computational  technique,  hereafter  called  the  embedded  well  technique, 
developed  to  compute  flow  rates  through  well  screens  is  described.  This  technique  is  applicable 
for  groundwater  flow  model  simulations  using  the  finite  element  method.  It  computes  for 
pumping  wells  and  relief  wells,  at  each  screen-associated  3D  mesh  node,  the  mass-conservative 
nodal  flow  (Cheng  et  al.  2010)  that  represents  the  flow  rate  through  the  well  screen  section 
associated  with  the  node.  The  technique  accounts  for  scenarios  that  include  or  ignore  the  head 
losses  across  the  well  screen  and  along  the  well.  The  technique  is  verified  with  a  simple  box 
model  that  contains  one  pumping  well  and  one  relief  well. 

METHODOLOGY:  To  solve  the  coupled  groundwater/well  flow  system  using  the  finite  element 
method,  it  is  convenient  to  represent  each  well  with  multiple  nodes  sitting  along  the  centerline  of 
the  well  (Konikow  et  al.  2009),  where  each  well  node  has  a  corresponding  subsurface  node 
sharing  the  same  coordinates,  as  shown  in  Figure  1.  In  the  figure,  the  well  is  represented  with 
nine  well  nodes,  where  the  top  five  nodes  are  associated  with  the  well  casing  and  the  bottom  four 
nodes  with  the  well  screen.  Each  well  node  is  associated  with  a  corresponding  3D  GW  node, 
where  GW  node  IDs  are  presented  in  red  color,  well  node  IDs  are  presented  in  blue  color,  and  the 
well  screen  is  highlighted  using  a  green  shade. 
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Figure  1 .  Example  of  a  multinode  well  and  its  associated  3D  GW  nodes. 

The  embedded  well  technique  computes  both  GW  heads  and  well  heads  at  model  mesh  nodes.  It 
solves  the  governing  equations  of  the  coupled  GW/well  system  that  include  the  3D  Richards 
equation  for  subsurface  flow  and  a  one-dimensional  (ID)  steady-state  equation  for  well  flow. 
The  Richards  equation  can  be  written  as 

^  +  V-VGW  =QGW  +Rgw  (1) 

dt 

where  0  is  moisture  content  [L3/L3],  1  is  time,  V  is  the  del  operator,  \JGW  is  Darcy  velocity 
[L/t],  QGW  is  the  source/sink  [L3/(L3t)]  of  the  GW  system  due  to  GW/well  interaction,  and  Raw 
represents  the  other  sources/sinks  [L3/(L3t)]  of  the  GW  system. 

With  the  Darcy’s  law,  VGW  =  — K6"  •  VHGW ,  Equation  1  can  be  further  written  as 

^_v.  Kgw  VHgw  =Qgw+Rgw  (2) 

dt 
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where  Kf;"  is  GW  hydraulic  conductivity  tensor  [L/t],  H GW  is  GW  total  head  [L], 

The  flow  pattern  within  a  well  can  be  complex  and  nonlinear  when  withdrawal  or  injection  is 
initiated.  To  resolve  the  temporal  variation  of  the  well  flow  pattern  accurately,  it  is  necessary  to 
solve  the  continuity  equation  and  the  momentum  equation  using  small  time-steps.  With  the 
assumption  that  the  well  flow  reaches  equilibrium  much  quicker  than  the  surrounding  local  GW 
flow,  a  ID  linear  steady-state  flow  equation  is  employed  to  represent  well  flow  in  the  coupled 
GW/well  system  so  that  much  greater  time-steps  can  be  used  for  computation.  The  ID  steady- 
state  well  equation  can  be  written  as 


t  well 


dV 


well 


Q"eU  +R 


dl 


well 


well 


(3) 


where  AweU  is  the  cross-sectional  area  of  the  well  [L2],  VweU  is  well  flow  velocity  [L/t],  lweU  is 
the  axis  along  well,  Q"eU  represents  the  source/sink  [L3/(Lt)]  of  the  well  system  due  to  GW/well 

11  Q 

interaction,  and  R"e  represents  the  other  source/sink  [L  /(Lt)]  of  the  well  system. 


Substituting  VweU 


-K 


well 


dHwe" 

dr" 


into  Equation  3  yields 


j^well 


d 


K 


well 


dHwe" 

dl"* 


dr" 


Qwell  _|_  j^well 


(4) 


where  KweU  represents  the  equivalent  hydraulic  conductivity  [L/t]  of  the  well  and  HweU  is  well 
total  head  [L], 

The  subsurface  flow  system  and  the  well  system  are  coupled  via  the  continuity  of  flux  that 
simply  states  that  water  leaving  the  GW  system  is  equal  to  water  entering  the  well  system,  and 
vice  versa,  at  any  well  locations  and  any  times,  as  shown  in  Equation  5. 

Q7"(t)  =  -Q7(t)  (5) 

where  Q'7"  (t)  is  the  net  flow  rate  [L3/t]  entering  the  well  at  a  well  location  K  and  time  t,  Q™  (t) 
is  the  net  flow  rate  [L3/t]  entering  the  GW  system  at  the  same  location  and  time. 

Finite  Element  Discretization.  Through  the  FE  discretization,  the  GW  equation  expressed 
with  Equation  2  can  be  approximated  (Lin  et  al.  1997)  as 

E-GWiGW  ,GW  ^ GW  ,  D GW  T  _ 

ai,KkK  ~bI  =Ql  +R,  I^NCW  I6) 


where  Ncw  is  the  total  number  of  3D  GW  nodes,  I  denotes  the  equation  associated  with  GW 
node  I,  K  denotes  a  GW  node  that  belongs  to  an  element  containing  node  I,  l\w  is  the  pressure 
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head  [L]  of  GW  node  K,  afj  and  bfw  are  coefficients  resulting  from  matrix  assembly  in  the 

finite  element  discretization,  Qfw  represents  the  nodal  source  [L3/t]  from  well,  and  Rfw 

represents  the  nodal  source  [L3/t]  via  the  other  means,  which  is  usually  zero  when  GW  node  I  is 
associated  with  the  well  screen. 

Likewise,  the  well  equation  expressed  with  Equation  4  can  be  approximated  as 

E  well  7  well  7  well  s~\well  ,  r% well  t  s-  at  Tri\ 

a.J,L  kL  ~bj  =Qj  +RJ  J  ^Nwell  (7) 

L 

where  Nwell  is  the  total  number  of  ID  well  nodes,  J  denotes  the  equation  associated  with  well 
node  J,  L  denotes  a  well  node  that  belongs  to  a  ID  element  around  node  J ,h™eU  is  the  pressure 
head  [L]  of  well  node  L,  a"re"  and  bJeU  are  coefficients  resulting  from  assembling  the  matrix  in 
the  finite  element  discretization,  QJ  represents  the  nodal  source  [L  /t]  from  subsurface,  and 
Ry  represents  the  nodal  source  [L  /t]  through  the  other  means  (e.g.,  withdrawal  from  a  water 
pump  placed  in  a  pumping  well  or  excess  groundwater  leaving  a  relief  well  from  top  of  the  well.) 

The  following  continuity  equations  thus  exist: 

qgw  _  _Qweii  _  q  if  we]]  no(ie  j  js  associated  with  the  casing; 

qgw  _  _Qjeii  ^  weji  node  j  is  associated  with  the  screen. 

where  GW  node  I  is  associated  with  well  node  J. 

Computing  Nodal  Flow.  When  the  Galerkin  FE  method  (Miller  et  al.  1998)  is  employed  to 
solve  the  Richards  equation,  the  residual  associated  with  each  node  can  be  computed  by 
substituting  the  computed  pressure  head  back  into  the  global  matrix  equation  (i.e.,  the  left-hand 
side  of  Equation  6)  that  is  constructed  from  the  FE  matrix  assembly  without  taking  into  account 
boundary  conditions  and  source/sink  terms  (Cheng  et  al.  2010).  If  the  node  is  an  internal  node, 
the  computed  residual  will  be  smaller  than  the  specified  residual  error  tolerance.  If  this  node  is  a 
boundary  node,  the  computed  residual  represents  the  net  flow  entering  or  leaving  the  GW 
computational  domain  through  the  boundary  face  area  associated  with  the  node.  If  this  node  is  a 

point  source  or  sink  in  the  domain,  then  the  computed  residual  is  the  injection  or  withdrawal  rate 

via  the  node.  Therefore,  the  computed  residual  of  a  screen-associated  GW  node  represents  the 
net  flow  from  GW  to  well  (when  the  residual  is  negative)  or  from  well  to  GW  (when  the  residual 
is  positive)  at  the  GW  node.  By  summing  up  the  flow  rates  of  all  GW  nodes  associated  with  the 
screen  of  a  well,  the  total  flow  rate  from  GW  to  well  is  thus  calculated,  as  illustrated  in  Figure  2 
where  Equation  8  is  employed  to  determine  QXwkeU  ,  Q2"kel1 ,  Q3"kel1 ,  Q4'keJI ,  and  Q5'kel1  at  each 

well  screen  node,  in  which  k  denotes  the  local  well  node  ID.  As  shown  in  Figure  2,  the  net  flow 
rate  from  GW  to  well  will  be  equal  to  the  rate  of  outflow  at  a  relief  well  or  the  specified  pumping 
rate  at  a  pumping  well,  as  a  steady  well  flow  is  assumed. 
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Flow  Zone  1 
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Q3 


-»Q2 


Q5 


Q4 
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Figure  2.  Concept  of  using  mass-conservative  GW  nodal  flow  to  compute  flow  rates  through  well 
screens  and  to  derive  the  outflows  of  relief  wells. 


Head  Loss  across  the  Well  Screen.  When  the  head  loss  across  the  well  screen  is  not 
negligible,  the  flow  rate  across  the  well  screen  can  be  computed  using  Equation  9: 


<2 


well 


qGW  =  j^well  _  screen 


GW 


■pwell 


H 


well 

J 


. screen 


(9) 


where  well  node  J  is  associated  with  the  well  screen,  GW  node  I  is  associated  with  well  node  J, 
Hfw  is  the  total  head  of  GW  node  I,  HjeU  is  the  total  head  of  well  node  J,  £jell-screen  |s  the 

equivalent  resistance  across  the  well  screen  associated  with  well  node  J,  AJell-screen  is  the 
equivalent  screen  area  associated  with  well  node  J  that  can  be  further  expressed  as 
2 KrweULwJell-screen  with  rweU  the  radius  of  the  well  and  Lw/u-screen  the  equivalent  length  of  the  well 

screen  associated  with  well  node  J,  provided  that  the  thickness  of  well  screen  packing  is 
relatively  small  when  compared  to  the  radius  of  the  well. 

Applying  the  Darcy’s  law  for  radial,  steady-state  flow  to  a  well  yields  Equation  10: 

Q  =  -2nKr -  (10) 

dr 
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where  H  is  total  head,  r  is  the  distance  in  the  radial  direction,  Q  is  the  flow  rate  in  the  radial 
direction,  and  K  is  the  hydraulic  conductivity  of  the  porous  medium  surrounding  the  well. 

By  assuming  a  radial  flow  through  the  well  screen,  Hfw  and  HJeU  in  Equation  9  can  be  related 
by  solving  Equation  10  as 


tjGW 

HI 


H 


well 


Qj 


well 


2  tzK 


welLf 


-In 


rweu  +0.5  d" 


well 


(11) 


where  K"eUscr“"  and  dwe“scmn  are  the  equivalent  hydraulic  conductivity  and  the  thickness, 
respectively,  of  well  screen  packing. 

Substituting  Equation  9  into  Equation  1 1  yields 


77  well  _  screen 

ej 


well 


•In 


rwell+0.5d 


well . 


well 


jrwell  _  screen 


(12) 


When  the  head  loss  across  the  well  screen  is  negligible,  the  resistance  is  very  small,  and 

j^well  _  screen 

J  n  becomes  very  large.  As  a  result,  the  algebraic  equations  of  the  coupled  system  may 

ET  ~screen 

become  ill-conditioned  and  difficult  to  solve.  To  overcome  this  difficulty,  the  total  head  at  a  well 
node  is  set  to  be  the  same  as  the  total  head  at  the  associated  GW  node  for  computation,  as  shown 
in  Equation  13: 


Hfw  =  HJe“  if  well  node  J  is  associated  with  a  screen.  (13) 


Head  Loss  along  the  Well.  When  the  head  loss  along  the  well  is  not  negligible,  Equation  4  is 
solved  to  compute  the  head  distribution  along  the  well.  To  solve  Equation  4,  the  well-known 
Hazen-Williams  equation  (i.e.,  Equation  14)  that  has  been  widely  used  to  estimate  the  head  loss 
in  steady  pipe  flow  can  be  used  to  estimate  KweU : 


hL=L 


Qpipe 


C„nit  '  EHW  A  - R 


.0.63 


1.852 


(14) 


where  hL  is  the  head  loss  in  steady  pipe  flow,  L  is  the  length  of  the  pipe,  Qpipe  is  the  steady 

volumetric  flow  rate,  A  is  the  cross-sectional  area  of  the  pipe,  R  is  the  hydraulic  radius  of  the 
pipe,  CIIW  is  the  Hazen-Williams  roughness  coefficient  (htty://www.  ensineerinstoolbox.  com/ 

hazen-williams-coefficients-d  798. html )  for  the  inside  wall  of  the  pipe,  and  Cunit  is  a  unit- 
specific  coefficient  in  the  Hazen-Williams  equation,  where  Cunit  is  1.318  for  U.S.  customary 
units  (foot  and  second)  and  0.849  for  SI  units  (meter  and  second). 

Equation  14  can  be  rewritten  as 
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t/  _  Qpipe  p0.63  y-7  ZJ  ^ unit  ^HW  ^ 

V  pipe  —  .  —  '^unit'^HW  ' R  '  V/l  —  ^  _  _  0.46 


.0.63 


■VH 


VH 


Therefore,  KweU  can  be  estimated  using  the  following  equation: 


D1 

Jewell  _  ^unit  J^HW  ‘  _ 


0.63 


VH 


0.46 


(15) 


(16) 


Because  the  Hazen-Williams  equation  is  applicable  when  the  flow  velocity  does  not  exceed  10 
ft/s,  one  may  use  this  velocity  value  to  calculate  K"'eU  when  there  are  not  field  data  available  for 
calibration.  By  assuming  laminar  flow,  one  may  also  use  Re  <  2,100  to  calculate  the  flow 

velocity  and  then  KweU ,  where  Re  is  the  Reynolds  number. 

Equation  4  would  require  small  time-steps  to  solve  when  KweU  is  large.  Because  the  head  loss 
along  the  well  becomes  negligible  when  KweU  is  sufficiently  large,  the  hydrostatic  condition,  as 
described  in  Equation  17,  can  be  employed  for  computation  so  that  reasonably  large  time-steps 
can  be  used.  The  entire  well  system  in  this  case  can  be  thought  of  as  a  super  mixer  within  which 
the  hydrostatic  condition  is  reached  immediately  under  any  disturbances. 

HweU=Hwell  (17) 

In  addition,  Equation  18  that  describes  mass  conservation  for  the  well  is  used  to  close  the 
computational  system: 


s-\  well 
^ outflow 


’  screen 

E 


Hfw  -H 


t  well  _  screen 


well 


-^w  ell  _  screen 


(18) 


where  HweU  is  the  hydrostatic  total  head  of  the  well,  and  Q''n^hnv  represents  the  rate  of  well 
outflow  to  balance  out  the  well  inflow  through  the  well  screen. 

Solving  the  Coupled  GW/Well  System.  Based  on  the  discussion  above  for  scenarios  with 
and  without  the  head  losses  across  the  well  screen  and  along  the  well,  various  scenarios  are 
discussed  as  follows.  To  help  discussion,  a  coupled  system  consisting  of  Now  groundwater  nodes 
and  one  well  is  used,  where  the  well  contains  Nweu  well  nodes,  which  includes  Nscreen  nodes 
associated  with  the  well  screen  and  Ncasing  nodes  associated  with  the  well  casing.  For  example,  in 
Figure  1,  Nwen  =  9,  Nscreen  =  4,  and  N casing  =  5.  For  the  coupled  systems  with  more  than  one  well, 
the  same  modeling  techniques  described  below  can  be  applied  directly. 

With  the  head  loss  across  the  screen :  When  the  head  loss  across  the  well  screen  is  accounted  for, 
the  GW  total  head  and  the  well  total  head  at  a  well  location  can  be  different.  In  this  case,  the  GW 
system  and  the  well  system  are  coupled  with  Equation  9. 
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If  the  head  loss  along  the  well  is  also  taken  into  account,  there  will  be  in  total  ( NGw  +  Nweii) 
unknowns  for  the  coupled  system,  which  includes  Ngw  GW  nodal  heads  and  Nweu  well  heads. 
The  equations  to  solve  include  Ngw  Equation  6  and  Nweu  Equation  7,  where  the  well  head  and  the 
associated  GW  head  at  a  screen  node  location  are  related  using  Equation  9.  Suppose  there  are 
1,200  GW  nodes  for  the  coupled  system  in  Figure  1,  the  Nweu  well  heads  are  treated  as  the  extra 
GW  unknowns  in  the  computational  system,  as  depicted  in  Figure  3  a. 

If  the  head  loss  along  the  well  is  negligible,  the  hydrostatic  condition  (i.e.,  Equation  17)  applies, 
and  there  is  actually  only  one  total  head  to  be  computed  for  the  entire  well.  In  this  case,  the 
coupled  system  has  in  total  (Ngw  +1)  unknowns,  as  illustrated  in  Figure  3b.  In  addition  to  the 
New  Equation  6,  Equation  1 8  is  employed  to  close  the  computational  system. 
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Figure  3.  Numbering  of  unknown  IDs  for  the  coupled 
GW/well  system  when  head  loss  across  the 
screen  is  present:  (a)  head  loss  along  the  well 
also  exists;  (b)  head  loss  along  the  well  is 
negligible. 


Without  the  head  loss  across  the  screen :  When  the  head  loss  across  the  well  screen  is  negligible, 
the  GW  total  head  is  equal  to  the  well  total  head  at  a  screen  location.  In  this  case,  the  GW  system 
and  the  well  system  are  coupled  with  Equation  13. 

If  the  head  loss  along  the  well  is  taken  into  account,  there  will  be  in  total  (New  +  Ncasing ) 
unknowns  for  the  coupled  system  because  the  well  head  and  the  associated  GW  head  at  each 
screen  node  location  are  identical.  In  this  case,  Equation  7  for  well  node  J  can  be  combined  with 
Equation  6  for  GW  node  I  to  form  the  following  equation: 
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E  GW  jGW 
ai  K  hK 


—  b 


GW 


—  bjdl  —  Rfw  +Rje" 

I  G  ^ GW  ,  j  G  NSI 


(19) 


where  GW  node  I  and  well  node  J  are  associated  with  the  same  well  location. 

As  a  result,  the  number  of  equations  will  be  the  same  as  the  number  of  unknowns.  Figure  4a 
shows  an  example  of  the  numbering  system  for  this  scenario. 


(a) 

GW  Well  Extra 
Node  ID  Node  ID  Unknown  ID 


145  A  1 


264 


357 


481 


622 


733 


858 


972  +  8 


1003 


1201 


1202 


1203 


1204 


1205 


Screened 

Section 


I 


(b) 

GW  Well 
Node  ID  Node  ID 

145 


264 


357 


481 


622 


733 


858 


972 


1003 


1 


Extra 
Unknown  ID 

1201 


Figure  4.  Numbering  of  unknown  IDs  for  the  coupled 
GW/well  system  when  head  loss  across  the 
screen  is  negligible:  (a)  head  loss  along  the  well 
exists;  (b)  head  loss  along  the  well  is  negligible. 


If  the  head  loss  along  the  well  is  negligible,  the  hydrostatic  condition  (Equation  17)  applies,  and 
there  is  only  one  total  head  needed  to  represent  the  entire  well  as  discussed  previously.  There  are 
thus  (Nscreen  +1)  unknowns  in  total  for  the  coupled  system,  as  illustrated  in  Figure  4b.  Because 

this  scenario  also  has  negligible  head  loss  across  the  well  screen,  all  GW  nodes  associated  with 
the  well  screen  are  set  to  have  a  total  head  identical  to  the  well  total  head.  The  entire  well  and  all 
the  GW  nodes  associated  with  the  well  screen  can  thus  be  thought  of  as  being  included  in  an 
imaginary  equalizer.  This  further  reduces  the  number  of  unknowns  to  (Nscreen  +1  —  Nscreen) .  As  a 

result,  the  net  rate  of  flow  entering  these  GW  nodes  is  also  the  net  rate  of  flow  from  the  GW 
system  to  the  well  system. 
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To  solve  the  coupled  system,  the  discretized  GW  equations  (i.e.,  Equation  6)  of  screen- 
associated  GW  nodes  (e.g.,  733,  858,  972,  and  1003  in  Figure  4b)  are  added  up  to  account  for  the 
interaction  between  the  imaginary  equalizer  and  the  surrounding  GW  system,  which  results  in 


E  E 


/c  yv 


GW 


E  a"  +  E 


GW 


IeN” 


leNz. 


(20) 


where  NfJeen  represents  the  set  of  GW  nodes  associated  with  the  well  screen. 

Equation  20  can  be  further  written  as  follows  using  mass  balance  within  the  well  at  a  steady  state 
(i.e.,  Y,  Qf"  +  Rj'e“  =0)  and  the  continuity  of  flux  through  the  GW-well  interface  (i.e., 

J^Nwell 

E  a“  =  -  E  er"  )■ 

TpA[GW  JEN 

1  screen  screen 


Jen\ 


■GW 

screen 


E  Eff  ~b°w  =  E  «r"  +  E 


GW 


(21) 


yew,„ 


/e# 


GW 

screen 


Equation  21  states  the  mass  balance  of  the  imaginary  equalizer  system.  To  close  the 
computational  system,  Equation  13  is  utilized  to  replace  the  N™een  Equation  6  of  well-associated 
GW  nodes. 


Boundary  conditions :  For  solving  the  3D  Richards  equation,  both  the  head-  and  the  flux-type 
boundary  conditions  can  be  applied.  For  instance,  the  embedded  well  technique  has  been 
incorporated  into  the  3D  subsurface  module  of  the  WASH123D  model  (Yeh  et  al.  2006), 
WASH3D  hereafter,  which  allows  the  modeler  to  apply  Dirichlet,  Cauchy,  Neumann,  and 
Variable  boundary  conditions  to  various  parts  of  the  domain  boundary  as  necessary.  As  for 
solving  the  ID  steady-state  well  flow  equation  which  is  embedded  in  the  3D  GW  computation,  a 
no-flow  boundary  condition  is  applied  to  both  ends  of  the  ID  well  domain,  as  the  pumping  rate 
of  a  pumping  well  and  the  outflow  rate  of  a  relief  well  are  treated  as  a  source/sink  term  to  their 
associated  well  systems  (i.e.,  RweU  in  Equation  4). 

Pump  location:  When  a  well  is  a  pumping  well,  the  location  of  the  water  pump  in  the  well  and 
the  pumping  schedule  (i.e.,  the  time-series  of  withdrawal/injection  rate)  are  given  for 
computation.  In  other  words,  RjeU  is  specified,  and  the  following  equation  exists: 


Y  Qfw  =  y  RT"  (22) 

JeNwdl 


As  mentioned  previously,  Qfw  cm  be  computed  by  substituting  the  computed  heads  into 
Equation  6,  with  Rfw  being  either  zero  or  specified. 

Relief  well:  Overflow  at  a  relief  well  may  or  may  not  occur,  depending  on  whether  excess  (i.e., 
high-pressure)  groundwater  exists  around  the  well  screen.  When  overflow  occurs,  the  total  head 
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at  well  top  is  controlled  at  user-specified  values  for  computation.  That  is,  water  is  assumed  to 
flow  freely  from  the  top  of  the  well,  and  the  pressure  head  at  the  overflow  location  is  assumed  to 
be  zero.  When  overflow  ceases,  the  well  water  level  is  below  that  specified  value  and  is 
computed  as  part  of  the  solution.  In  this  case,  there  is  no  flow  calculated  at  well  nodes  above  the 
well  water  level,  and  the  heads  at  those  well  nodes  will  be  undefined  and  excluded  from 
computation  because  they  are  excluded  from  the  computation. 

Figure  5  depicts  the  flow  chart  of  the  computational  algorithm  for  handling  relief  wells  in  the 
embedded  well  technique.  As  shown  in  the  figure,  for  any  single-step  computation,  which  is  either 
a  steady-state  computation  or  a  transient  simulation  over  one  time-step,  the  user-specified  total 
head  values  applied  to  relief  wells  are  first  examined.  For  relief  wells  where  the  specified  total 
head  values  are  greater  than  the  given  reference  elevations,  the  specified  total  head  values  apply  to 
the  top  of  those  relief  wells  as  Dirichlet  boundary  conditions,  and  they  are  flagged  as  active  relief 
wells.  Otherwise,  no  boundary  condition  is  applied,  and  the  relief  well  is  flagged  as  inactive  even 
though  flow  is  free  to  move  in  and  out  of  the  well.  After  the  coupled  GW/well  system  is  solved  for 
the  first  time,  the  computed  head  value  (i.e.,  the  initial  solution  in  Figure  5)  at  the  top  of  each  relief 
well  flagged  inactive  is  checked.  If  the  computed  total  head  is  greater  than  the  given  reference 
elevation  of  an  inactive  relief  well,  the  relief  well  is  flagged  active.  Then  the  given  reference 
elevation  is  applied  to  the  top  of  the  well  as  a  Dirichlet  boundary  condition  in  solving  the  coupled 
system  for  the  second  time.  The  coupled  system  needs  to  be  solved  the  second  time  only  when  any 
relief  wells  switch  from  inactive  to  active  based  on  the  initial  solution. 


Figure  5.  Flow  chart  to  solve  the  coupled  GW/well  system  when  relief  wells 
exist. 
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Input  Data.  The  WASH3D  model  uses  card  input  to  read  data  needed  for  user-specified 
simulation.  To  account  for  embedded  wells,  RW1  and  RW3  cards  have  been  added  to  the  input 
file  that  defines  embedded  relief  and  pumping  wells,  respectively,  in  subsurface  flow 
simulations.  An  RWO  card  defines  a  maximum  saturated  hydraulic  conductivity  value  allowed 
for  computing  ID  steady  flow  and  a  minimum  well  screen  resistance  coefficient  value  allowed 
for  computing  head  loss  across  the  well  screen.  The  contents  of  these  cards  are  defined  in  Tables 
1  through  3. 


Table  1.  Contents  of  RWO  card  for  embedded  wells:  2  parameters  are  included. 

Field 

Type 

Parameter  Name 

Parameter  Description 

1 

Character 

RWO 

Card  identifier  for  parameters  for  all  embedded  wells. 

2 

Positive  Real 

AK_RWF_MAX 

Maximum  saturated  hydraulic  conductivity  value  allowed 
for  computing  1 D  steady  flow. 

3 

Positive  Real 

RRWFMIN 

Minimum  well  screen  resistance  coefficient  value  allowed 
for  computing  head  loss  across  the  well  screen. 

Table  2.  Contents  of  RW1  card  for  embedded  relief  wells:  9  parameters  are  included. 

Field 

Type 

Parameter  Name 

Parameter  Description 

1 

Character 

RW1 

Card  identifier  for  embedded  relief  well. 

2 

Positive  Integer 

IDRWF1  (1) 

ID  of  the  embedded  well  corresponding  to  the  l-th 
sequentially-listed  embedded  well  that  is  a  relief  well. 

3 

Positive  Integer 

NRWNF(I) 

Total  number  of  well  nodes  associated  with  the  l-th 
sequentially  listed  embedded  well. 

4 

Positive  Integer 

NRWSNF(I) 

Number  of  well  screen  nodes  associated  with  the  l-th 
sequentially  listed  embedded  well  feature  that  is  a  relief 
well. 

5 

Positive  Integer 

IRWTYPF(I) 

ID  of  the  x-y  series  to  describe  the  time-dependent  head 
[L]  controlled  at  top  of  the  well  for  the  l-th  sequentially 
listed  embedded  well. 

6 

Integer  (0  or  1) 

IDBH(I) 

Index  of  the  profile  type  for  the  l-th  sequentially  listed 
embedded  well:  IDBH(I)  =  0,  if  it  is  a  pressure  head  profile; 
IDBH(I)  =  1 ,  if  it  is  a  total  head  profile. 

7 

Positive  Real 

Z_RWF(I) 

Reference  elevation  of  the  l-th  sequentially  listed 
embedded  well,  which  is  used  to  indicate  whether  overflow 
occurs  at  a  relief  well. 

8 

Positive  Real 

AK_RWF(I) 

Saturated  hydraulic  conductivity  of  the  l-th  sequentially 
listed  embedded  well. 

9 

Positive  Real 

D_RWF(I) 

Well  diameter  of  the  l-th  sequentially  listed  embedded 
well. 

10 

Positive  Real 

R_RWF(I) 

Well  screen  resistance  coefficient  associated  with  the  l-th 
sequentially  listed  embedded  well. 
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Table  3.  Contents  of  RW3  card  for  embedded  pumping  wells:  9  parameters  are  included. 

Type 

Parameter  Name 

Parameter  Description 

1 

Character 

RW3 

Card  identifier  for  embedded  pumping  well. 

2 

Positive  Integer 

IDRWF1  (1) 

ID  of  the  embedded  well  corresponding  to  the  I-th 
sequentially  listed  embedded  well  that  is  a  relief  well. 

3 

Positive  Integer 

NRWNF(I) 

Total  number  of  well  nodes  associated  with  the  I-th 
sequentially  listed  embedded  well. 

4 

Positive  Integer 

NRWSNF(I) 

Number  of  well  screen  nodes  associated  with  the  I-th 
sequentially  listed  embedded  well  that  is  a  relief  well. 

5 

Positive  Integer 

IRWTYPF(I) 

ID  of  the  x-y  series  to  describe  the  time-dependent 
pumping  rate  [L3/t]  for  the  I-th  sequentially  listed 
embedded  well. 

6 

Positive  Integer 

IDRWF3(I) 

Location  of  the  pump  placed  in  the  I-th  sequentially  listed 
embedded  well  that  is  a  PUMPING  well.  For  example,  if 
IDRWF3(3,I)  =  5,  this  pump  is  placed  at  the  5-th  well  node 
from  top. 

7 

Positive  Real 

Z_RWF(I) 

Reference  elevation  of  the  I-th  sequentially  listed 
embedded  well,  which  is  used  to  indicate  whether 
overflow  occurs  at  a  relief  well. 

8 

Positive  Real 

AK_RWF(I) 

Saturated  hydraulic  conductivity  of  the  I-th  sequentially 
listed  embedded  well. 

9 

Positive  Real 

D_RWF(I) 

Well  diameter  of  the  I-th  sequentially  listed  embedded 
well. 

10 

Positive  Real 

R_RWF(I) 

Well  screen  resistance  coefficient  associated  with  the  I-th 
sequentially  listed  embedded  well. 

When  a  saturated  hydraulic  conductivity  value  (i.e.,  AK_RWF(I))  given  in  the  RW1  or  the  RW3 
card  is  greater  than  the  maximum  saturated  conductivity  allowed  (i.e.,  AK_RWF_MAX)  given 
in  the  RWO  card,  head  loss  along  the  well  is  assumed  negligible,  and  the  hydrostatic  condition  is 
applied  to  the  well  for  computation,  whether  it  is  a  relief  well  or  a  pumping  well.  Likewise,  when 
a  well  screen  resistance  coefficient  value  (i.e.,  R_RWF(I))  is  smaller  than  the  given  minimum 
resistance  coefficient  allowed  (i.e.,  R_RWF_MIN),  head  loss  across  the  well  screen  is  assumed 
negligible,  and  the  total  head  at  a  well  screen  node  is  considered  identical  to  that  of  its 
corresponding  GW  node. 

To  identify  the  groundwater  nodes  associated  with  each  well  and  its  screen  section(s),  two  more 
lines  are  needed  after  each  RW1  or  RW3  card,  and  each  line  contains  NRWNF(I)  record,  where 
NRWNF(I)  is  the  total  number  of  well  nodes  associated  with  the  I-th  sequentially  listed 
embedded  well.  The  first  line  lists  the  IDs  of  3D  GW  global  nodes  associated  with  the  I-th 
sequentially  listed  embedded  well.  These  nodes  are  listed  in  order  from  top  down  along  the  well 
(i.e.,  the  first  node  is  the  top  well  node,  and  the  last  node  is  the  bottom  well  node).  The  second 
line  defines  nodes  associated  with  the  well  screen  and  well  casing  of  the  embedded  well.  A 
screen  node  is  indicated  with  an  integer  value  of  1,  while  a  casing  node  is  given  a  value  of  0. 
These  values  are  also  assigned  to  the  well  nodes  from  top  down  along  the  well. 

Below  is  an  example  with  both  RW1  and  RW3  cards  for  two  embedded  wells  in  a  3D  GW 
model,  where  RW1  and  RW3  cards  define  a  relief  well  and  a  pumping  well,  respectively.  As 
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given  below,  both  wells  are  discretized  using  nine  nodes,  where  the  relief  well  has  two  screen 
sections  with  three  nodes  included  in  each  section,  and  the  pumping  well  has  only  one  screen 
section  that  covers  the  bottom  three  nodes.  Both  wells  account  for  head  losses  across  the  well 
screen  and  along  the  well  based  on  the  well  conductivity  values  and  the  well  screen  resistance 
coefficient  values  given  in  fields  8  and  10,  respectively. 

RWO  9999.999  0.00001 
RW1  1  9  6  10  1  10.0  2000.0  0.5  0.001 

173  591  800  1009  1218  1427  1636  1845  2054 

001110111 

RW3  2  9  3  12  9  10.0  5000.0  3.0  0.001 

67  485  694  903  1112  1321  1530  1739  1948 
000000111 

Preprocessing.  A  preprocessor  was  developed  to  achieve  memory  allocation  for  the  specified 
simulation  as  well  as  to  generate  node  connectivity  information  for  the  FE  computation.  Based  on 
the  data  given  with  RWO,  RW1,  and  RW3  cards  described  previously,  the  preprocessor  determines 
the  extra  unknowns  needed  to  account  for  each  embedded  well,  which  combines  with  the  element 
indices  information  provided  in  the  geometry  file  to  generate  node  connectivity  information.  Table 
4  lists  the  numbers  of  extra  unknowns  associated  with  various  head  loss  scenarios  for  an  embedded 
well,  where  the  well  is  assumed  a  relief  well  in  Figure  1.  Figure  6  depicts  the  3D  GW  nodes  and 
extra  unknowns  associated  with  the  embedded  relief  well  for  the  four  scenarios  listed  in  Table  4, 
provided  that  there  are  1,200  nodes  in  the  3D  computational  mesh. 


Table  4.  Number  of  extra  3D  unknowns  needed  for  an  embedded  relief  well  at  various 
scenarios  regarding  head  losses. 

Scenario 

ID 

RW1  Card  Input 

Head 

Loss 

across 

Screen 

Head 

Loss 

along 

Well 

Number  of 
Extra 

Unknowns 

1 

RW1  1  9  4  10  1  1 0.0  2000.0  0.5  0.01 

145  264  357  481  622  733  858  972  1003 

000001 1 1 1 

Yes 

Yes 

9 

2 

RW1  1  9  6  10  1  1 0.0  20000.0  0.5  0.01 

145  264  357  481  622  733  858  972  1003 

000001 1 1 1 

Yes 

No 

1 

3 

RW1  1  9  6  10  1  1 0.0  2000.0  0.5  0.000005 

145  264  357  481  622  733  858  972  1003 

000001 1 1 1 

No 

Yes 

5 

4 

RW1  1  9  6  10  1  1 0.0  20000.0  0.5  0.000005 

145  264  357  481  622  733  858  972  1003 

000001 1 1 1 

No 

No 

1 
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Figure  6.  An  example  demonstrating  3D  GW  nodes,  1 D  well  nodes,  and  extra  unknowns  associated 
with  an  embedded  relief  well  at  various  scenarios  regarding  head  losses. 


EXAMPLE:  A  simple  box  model  (Figure  7)  was  constructed  to  verify  the  proposed  relief  well 
algorithm. 

Model  Setup.  As  shown  in  Figure  7,  the  computational  domain  had  a  volume  of 
100  x  100  x  10  ft3  (x  =  0  to  100  ft,  y  =  0  to  100  ft,  and  z  =  0  to  10  ft)  and  was  discretized  with 
2,299  nodes  and  3,280  elements.  The  3D  mesh  was  composed  of  both  hexahedral  and  triangular 
prism  elements.  This  model  had  the  same  Dirichlet  boundary  conditions  assigned  on  the  right 
and  left  boundary  surfaces  to  mimic  the  water  levels  of  the  nearby  lake/river  and  variable 
boundary  conditions  on  the  top  boundary  surface  to  account  for  rainfall  contribution.  This  box 
model  included  two  embedded  wells:  one  relief  well  and  one  pumping  well  (Figure  7). 

A  transient  simulation  of  10  days  was  conducted  with  the  initial  total  head  set  to  10  ft  for  the 
entire  domain.  The  rainfall  rate  was  zero  during  the  first  5  days  and  2.54  x  10"2  ft/day  during  the 
last  5  days.  The  water  depth  on  top  of  the  relief  well  was  set  to  0  m  during  the  first  5  days  when 
excess  groundwater  came  out  of  the  well,  and  a  water  depth  of  0.1  ft  was  controlled  during  the 
last  5  days.  The  variation  of  total  heads  applied  to  the  left  and  right  boundary  faces  and  the 
pumping  rate  specified  at  the  pumping  well  are  depicted  in  Figure  8.  The  domain  contained  three 
geologic  material  types  with  saturated  hydraulic  conductivity  values  set  to  30,  0.0001,  and 
1,000  ft/day,  respectively.  The  modified  compressibility  was  0.001  for  all  the  three  material 
types.  A  simple  set  of  soil  curves  for  both  material  types  was  employed  as  shown  in  Figure  9. 
The  time  interval  used  for  the  transient  simulation  was  0.5  day  throughout  the  simulation. 
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Figure  7.  Plan  view  (left)  and  oblique  view  (right)  of  the  box  model,  where  z-magnification  =  10. 


Figure  8.  Head  boundary  conditions  (left)  and  specified  pumping  rate  (right)  used  for  the  box  model. 
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Figure  9.  Soil  curves  used  for  the  box  model. 
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Results  and  Discussion.  Four  scenarios,  associated  with  different  combinations  of  with  and 
without  head  losses  along  wells  and  across  well  screens,  were  simulated.  Table  5  presents  the 
RW1  and  RW3  card  input  data  for  these  scenarios,  where  the  saturated  hydraulic  conductivity 
for  well  flow  value  (highlighted  in  red  color)  and  the  well  screen  resistance  coefficient  value 
(highlighted  in  blue  color)  were  compared  with  the  given  AK_RWF_MAX  (set  to  9999.999)  and 
R_RWF_MIN  (set  to  0.00001)  values,  respectively,  to  determine  whether  head  loss  is  negligible. 


Table  5.  RW1  and  RW3  card  input  data  for  scenarios  regarding  head  losses. 

Scenario 

ID 

RW1/RW3  Card  Input 

Head  Loss 
Across 

Screen 

Head  Loss 
Along  Well 

1 

RW1  1  9  6  10  1  10.0  2000.0  0.5  0.001 

173  591  800  1009  1218  1427  1636  1845  2054 
001110111 

RW3  2  9  3  12  9  10.0  5000.0  3.0  0.001 

67  485  694  903  1 1 1 2  1 321  1 530  1 739  1 948 

0000001 1 1 

Yes 

Yes 

2 

RW1  1  9  6  10  1  10.0  20000.0  0.5  0.001 

173  591  800  1009  1218  1427  1636  1845  2054 
001110111 

RW3  2  9  3  12  9  10.0  50000.0  3.0  0.001 

67  485  694  903  1 1 1 2  1 321  1 530  1 739  1 948 

0000001 1 1 

Yes 

No 

3 

RW1  1  9  6  10  1  10.0  2000.0  0.5  0.000001 

173  591  800  1009  1218  1427  1636  1845  2054 
001110111 

RW3  2  9  3  12  9  10.0  5000.0  3.0  0.000001 

67  485  694  903  1112  1321  1530  1739  1948 

0000001 1 1 

No 

Yes 

4 

RW1  1  9  6  10  1  10.0  20000.0  0.5  0.000001 

173  591  800  1009  1218  1427  1636  1845  2054 
001110111 

RW3  2  9  3  12  9  10.0  50000.0  3.0  0.000001 

67  485  694  903  1 1 1 2  1 321  1 530  1 739  1 948 

0000001 1 1 

No 

No 

Tables  6  and  7  compare  total  head  distribution  at  well  nodes  and  their  associated  GW  nodes  at 
time  =  5  days  and  10  days,  respectively.  Figure  10  provides  a  visual  comparison  of  GW  total 
head  distribution  and  flow  directions  around  the  two  wells  at  both  time  =  5  days  and  10  days. 
Table  8  compares  the  nodal  flow  rates,  from  GW  to  well,  at  well  nodes  at  time  =  5  days  and  10 
days,  where  the  total  well  flow  rate  is  also  listed. 
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Table  6.  Computed  well  and  the  associated  GW  total  heads  at  time  =  5  days. 

Well 

Node 

ID 

Attribute1 

Scenario  1z 

Scenario  22 

Scenario  32 

Scenario  42 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

RW1-1 

CN 

N/A3 

9.40 

N/A 

9.69 

N/A 

9.42 

N/A 

9.76 

RW1-2 

CN 

9.43 

9.41 

9.81 

9.70 

9.45 

9.43 

9.83 

9.77 

RW1-3 

SN 

9.43 

9.43 

9.81 

9.75 

9.45 

9.45 

9.83 

9.83 

RW1-4 

SN 

9.46 

9.46 

9.81 

9.78 

9.48 

9.48 

9.83 

9.83 

RW1-5 

SN 

9.52 

9.50 

9.81 

9.79 

9.53 

9.53 

9.83 

9.83 

RW1-6 

CN 

9.66 

9.66 

9.81 

9.62 

9.68 

9.65 

9.83 

9.62 

RW1-7 

SN 

9.81 

9.84 

9.81 

9.83 

9.84 

9.84 

9.83 

9.83 

RW1-8 

SN 

9.84 

9.84 

9.81 

9.83 

9.84 

9.84 

9.83 

9.83 

RW1-9 

SN 

9.84 

9.84 

9.81 

9.83 

9.84 

9.84 

9.83 

9.83 

RW3-1 

CN 

15.92 

9.17 

16.48 

9.18 

14.49 

9.17 

14.68 

9.18 

RW3-2 

CN 

15.92 

9.17 

16.48 

9.18 

14.49 

9.17 

14.68 

9.18 

RW3-3 

CN 

15.92 

9.17 

16.48 

9.18 

14.49 

9.17 

14.68 

9.18 

RW3-4 

CN 

15.92 

9.17 

16.48 

9.18 

14.49 

9.17 

14.68 

9.18 

RW3-5 

CN 

15.92 

9.17 

16.48 

9.18 

14.49 

9.17 

14.68 

9.18 

RW3-6 

CN 

15.92 

10.30 

16.48 

10.30 

14.49 

10.24 

14.68 

10.26 

RW3-7 

SN 

15.92 

15.02 

16.48 

14.99 

14.49 

14.49 

14.68 

14.68 

RW3-8 

SN 

16.16 

14.85 

16.48 

14.89 

14.56 

14.56 

14.68 

14.68 

RW3-9 

SN 

16.75 

14.52 

16.48 

14.50 

14.78 

14.78 

14.68 

14.68 

1  CN  =  casing  node;  SN  =  screen  node. 

2  Scenario  1 :  head  loss  along  the  well  and  across  the  screen; 


Scenario  2:  head  loss  across  the  screen  only; 

Scenario  3:  head  loss  along  the  well  only; 

Scenario  4:  negligible  head  losses. 

3  N/A  denotes  the  well  water  level  is  lower  than  the  elevation  of  the  well  node. 


As  can  be  seen  from  Tables  6  and  7,  the  water  level  at  a  well  node  (i.e.,  well  head)  was  the  same 
as  the  total  head  of  the  associated  GW  node  when  the  well  node  was  located  within  a  screen 
section  with  negligible  head  loss  across  the  screen  (i.e.,  Scenarios  3  and  4).  The  well  head  was 
different  from  the  GW  total  head  of  the  associated  GW  node  when  either  the  well  node  was  in  a 
casing  section  or  head  loss  across  the  screen  was  not  negligible  (i.e.,  Scenarios  1  and  2). 

At  time  =10  days,  the  high  head  applied  to  the  boundary  associated  with  layer  3  resulted  in 
overflow  at  the  relief  well,  where  a  controlled  water  level  at  10.1  ft  was  set  to  the  top  of  the  well 
(i.e.,  well  head  at  RW1-1,  Table  7).  The  differences  of  GW  total  head  distribution  between  the 
two  wells  were  more  obvious  among  scenarios,  when  compared  with  those  at  time  =  5  days 
(Figure  10).  In  Table  8,  the  outflow  from  the  relief  well  also  varied  drastically  among  these  four 
scenarios,  where  the  application  of  hydrostatic  condition  in  the  well  (i.e.,  negligible  head  loss 
along  well,  Scenarios  2  and  4)  would  yield  much  more  water  drawn  from  GW  to  the  well,  when 
compared  with  the  scenarios  accounting  for  head  loss  along  the  well  (i.e.,  Scenarios  1  and  3). 
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Table  7.  Computed  well  and  the  associated  GW  total  heads  at  time  =  10  days. 

Well 

Node 

ID 

Attribute1 

Scenario  1z 

Scenario  22 

Scenario  32 

Scenario  42 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

Well 
Head,  ft 

GW 

Head,  ft 

RW1-1 

CN 

10.10 

11.16 

10.10 

9.91 

10.10 

11.30 

10.10 

9.98 

RW1-2 

CN 

10.75 

11.23 

10.10 

9.93 

10.83 

11.36 

10.10 

10.01 

RW1-3 

SN 

11.39 

11.42 

10.10 

10.01 

11.56 

11.56 

10.10 

10.10 

RW1-4 

SN 

11.93 

11.87 

10.10 

10.05 

12.08 

12.08 

10.10 

10.10 

RW1-5 

SN 

12.71 

12.43 

10.10 

10.08 

12.77 

12.77 

10.10 

10.10 

RW1-6 

CN 

14.65 

10.91 

10.10 

10.36 

14.87 

10.95 

10.10 

10.10 

RW1-7 

SN 

16.58 

16.98 

10.10 

13.62 

16.96 

16.96 

10.10 

10.10 

RW1-8 

SN 

16.93 

17.00 

10.10 

13.87 

16.99 

16.99 

10.10 

10.10 

RW1-9 

SN 

16.99 

17.02 

10.10 

14.16 

17.01 

17.01 

10.10 

10.10 

RW3-1 

CN 

11.05 

9.09 

10.04 

8.96 

12.49 

9.11 

11.26 

8.98 

RW3-2 

CN 

11.05 

9.09 

10.04 

8.95 

12.49 

9.11 

11.26 

8.96 

RW3-3 

CN 

11.05 

9.09 

10.04 

8.95 

12.49 

9.11 

11.26 

8.96 

RW3-4 

CN 

11.05 

9.09 

10.04 

8.95 

12.49 

9.11 

11.26 

8.96 

RW3-5 

CN 

11.05 

9.09 

10.04 

8.95 

12.49 

9.11 

11.26 

8.96 

RW3-6 

CN 

11.05 

11.23 

10.04 

11.18 

12.49 

11.15 

11.26 

11.06 

RW3-7 

SN 

11.05 

11.95 

10.04 

11.53 

12.49 

12.49 

11.26 

11.26 

RW3-8 

SN 

10.81 

12.13 

10.04 

11.63 

12.42 

12.42 

11.26 

11.26 

RW3-9 

SN 

10.22 

12.46 

10.04 

12.02 

12.19 

12.19 

11.26 

11.26 

1  CN  =  casing  node;  SN  =  screen  node. 

2  Scenario  1 :  head  loss  along  the  well  and  across  the  screen; 
Scenario  2:  head  loss  across  the  screen  only; 

Scenario  3:  head  loss  along  the  well  only; 

Scenario  4:  negligible  head  losses. 


To  verify  the  computed  flow  rates  via  the  relief  well  were  calculated  correctly,  a  second  set  of 
WASH3D  simulations  was  conducted,  where  the  GW  nodes  associated  with  well  screens  were 
treated  as  point  sources/sinks,  and  their  respective  flow  rates  computed  previously  (e.g.,  those 
listed  in  Table  8)  were  employed  as  the  injection/withdrawal  rates.  Although  not  shown  here,  the 
computed  GW  total  heads  were  the  same  as  those  computed  previously  (e.g.,  those  listed  in 
Tables  6  and  7),  which  indicates  that  the  relief  well  flow  rates  were  computed  accurately. 

SUMMARY:  This  technical  note  describes  the  details  of  an  embedded  well  technique 
developed,  based  on  the  finite  element  method,  to  compute  mass-conservative  nodal  flow  rates  at 
screen  nodes  of  both  relief  and  pumping  wells.  The  technique  solves  the  governing  equations  of 
the  coupled  system  of  GW  and  embedded  wells  to  compute  the  nodal  flow  at  each  well  node 
location,  where  no  presumed  nodal  flow  distribution  over  well  screen  is  necessary.  As  a  result, 
the  GW  model  can  be  better  calibrated  and  the  outflow  of  a  relief  well  can  be  accurately 
estimated.  This  numerical  technique  accounts  for  scenarios  with  and  without  well  inefficiency 
(i.e.,  head  losses  along  the  well  and  across  the  screen).  This  technique  and  an  algorithm  to  handle 
the  computation  associated  with  relief  wells  have  been  incorporated  into  the  WASH3D  model 
and  can  be  incorporated  into  other  FEM-based  groundwater  models.  A  box  model  was  employed 
to  verify  the  implementation  of  embedded  technique  in  WASH3D. 
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Time  =  5  days 


Time  =  10  days 


Scenario  4: 
Negligible  Losses 


Total  Head 
”  18.0 
H  16.0 

14.0 
12.0 

110.0 
8.0 


Figure  10.  GW  total  head  distribution  (color  code)  and  flow  directions  (arrows)  around  wells  at  time 
5  days  (left  images)  and  time  =  10  days  (right  images),  where  z-magnification  =  5. 
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Table  8.  Computed  flow  rates  of  wells  at  time  =  5  da’ 

/s  and  time  =  10  days. 

Well  Node 
ID 

Scenario  I1 

Scenario  21 

Scenario  31 

Scenario  41 

Nodal  Flow2,  ft3/d 

Nodal  F 

ow,  ft3/d 

Nodal  F 

ow,  ft3/d 

Nodal  FI 

ow,  ft3/d 

5  days 

1 0  days 

5  days 

10  days 

5  days 

1 0  days 

5  days 

10  days 

RW1-1 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW1-2 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW1-3 

-10.82 

43.60 

-44.08 

-72.73 

-11.47 

82.60 

-71.61 

-115.64 

RW1-4 

-11.04 

-97.26 

-44.03 

-71.95 

-7.82 

-69.60 

-38.58 

-59.18 

RW1-5 

-35.73 

-451.66 

-28.81 

-46.89 

-42.39 

-548.70 

-19.26 

-29.56 

RW1-6 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW1-7 

47.24 

622.98 

55.71 

8295.55 

60.77 

809.81 

17.95 

5406.38 

RW1-8 

8.59 

113.21 

39.82 

5928.93 

0.37 

4.95 

36.02 

10837.20 

RW1-9 

1.76 

23.19 

21.39 

3185.09 

0.54 

7.18 

75.48 

22805.20 

Total  Well 
Flow3  (ft3/d) 

0.00 

254.06 

0.00 

17218.00 

0.00 

286.24 

0.00 

38844.30 

Well  Node 
ID 

Nodal  F 

ow,  ft3/d 

Nodal  F 

ow,  ft3/d 

Nodal  F 

ow,  ft3/d 

Nodal  F 

ow,  ft3/d 

5  days 

1 0  days 

5  days 

10  days 

5  days 

1 0  days 

5  days 

10  days 

RW3-1 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-2 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-3 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-4 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-5 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-6 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

RW3-7 

-8478.32 

8478.32 

-7032.01 

7032.01 

-2558.57 

2558.57 

-4361.57 

4361 .57 

RW3-8 

-12400.00 

12400.00 

-15043.20 

15043.20 

-5339.13 

5339.13 

-8743.83 

8743.83 

RW3-9 

-10537.58 

10537.58 

-9340.69 

9340.69 

-23518.20 

23518.20 

-18310.40 

18310.50 

Total  Well 
Flow3  (ft3/d) 

-31415.90 

31415.90 

-31415.90 

31415.90 

-31415.90 

31415.90 

-31415.90 

31415.90 

1  Scenario  1 :  head  loss  a 

ong  the  we 

1  and  across  the  screen; 

Scenario  2:  head  loss  across  the  screen  only; 

Scenario  3:  head  loss  along  the  well  only; 

Scenario  4:  negligible  head  losses. 

2  A  positive  nodal  flow  rate  denotes  water  entering  well  from  GW  at  the  well  node,  and  a  negative  value 
indicates  water  is  from  well  to  GW  at  the  well  node. 

3  For  a  relief  well,  a  positive  total  well  flow  represents  the  rate  of  outflow  at  the  well  top;  for  a  pumping 
well,  a  positive  total  flow  represents  its  withdrawal  rate,  and  a  negative  total  flow  represents  its  injection 
rate. 


ADDITIONAL  INFORMATION:  This  work  was  sponsored  by  the  Flood  &  Coastal  Storm 
Damage  Reduction  Program.  Dr.  Stacy  E.  Howington  and  Dr.  Matthew  W.  Farthing  of  the 
U.S.  Army  Engineer  Research  and  Development  Center  (ERDC),  Coastal  and  Hydraulics 
Laboratory  (CHL)  provided  valuable  inputs.  For  additional  information,  contact  Hwai-Ping 
(Pearce)  Cheng,  ERDC-CHL,  3909  Halls  Ferry  Road,  Vicksburg,  MS  39180,  at  601-634-3699  or 
e-mail:  Hwai-Pins.  Chens @  usace.  army.mil.  This  CHETN  should  be  cited  as  follows: 
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